
rm(list=ls(all=TRUE))
library(sf)
library(tidyverse)
library(cowplot)
library(readstata13)
library(openxlsx)
library(RColorBrewer)
library(survival)
library(plyr)

setwd(dirname(rstudioapi::getSourceEditorContext()$path))


####Figure 1#######################################################

rm(list=ls(all=TRUE))

county<-read_sf("./map/county.shp")
province<-read_sf("./map/province.shp")
line<-read_sf("./map/HSRline.shp")
China <- st_read("./map/China.json")

mapdata<-read.dta13("1data_replicate.dta")

data<-mapdata %>% filter(year==2004,general_living>0) %>%
  select(provinceid,prefectureid,countyid,cname,general_living) %>%
  rename(PAC=countyid)

data$PAC[data$PAC==130182]<-130109
data$PAC[data$PAC==350822]<-350803
data$PAC[data$PAC==360721]<-360704
data$PAC[data$PAC==370282]<-370215
data$PAC[data$PAC==371081]<-371003
data$PAC[data$PAC==371523]<-371503
data$PAC[data$PAC==430122]<-430112
data$PAC[data$PAC==511721]<-511703
data$PAC[data$PAC==610623]<-610681
data$PAC[data$PAC==610721]<-610703

data<-full_join(county,data)

table(data$general_living)
data$ncut <- cut(data$general_living, breaks = c(0,1,2,3,4,7,13),
                 labels=c('1','2','3','4','[5,7]','[8,13]'))

colourPalette <- brewer.pal(7,'Greens')[2:7]

ChinaMap<-ggplot() +
  geom_sf(data = China, colour = "grey",lwd =0.4,alpha=0.5)+
  geom_sf(data = data, aes(fill=ncut),colour = "white",alpha=0.8)+
  geom_sf(data = province,fill=NA,colour = "grey",lwd =0.4,alpha=0.5)+
  geom_sf(data = line,colour="black",linewidth=1)+
  geom_sf(data = line,colour="white",linetype = "dashed",linewidth=0.5)+
  coord_sf(ylim = c(-2387082,1654989),xlim = c(-2387082,2554989),crs="+proj=laea +lat_0=40 +lon_0=104")+
  scale_fill_manual(name="",values=colourPalette,na.value="white",na.translate=F)+
  theme_minimal()+
  theme(legend.position = c(1, 0.55),
        panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

SouthChinaSea <- ggplot() +
  geom_sf(data = China,  colour = "grey",lwd =0.4,alpha=0.5,fill=NA)+
  coord_sf(xlim = c(117131.4,2115095),ylim = c(-4028017,-1877844),
           crs="+proj=laea +lat_0=40 +lon_0=104")+
  theme(aspect.ratio = 1.25,axis.text = element_blank(),axis.ticks = element_blank(),
        axis.title = element_blank(),panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour="black",fill = NA,linewidth=1.5),
        plot.margin=unit(c(0,0,0,0),"mm"))

ppi <- 600
png("figure1.png", width=10*ppi, height=6*ppi, res=ppi)
ggdraw() +
  draw_plot(ChinaMap) +
  draw_plot(SouthChinaSea, x = 0.8, y = 0.00, width = 0.19, height = 0.35)
dev.off()


####Figure 3#######################################################

rm(list=ls(all=TRUE))
data<-read.xlsx("2data_figure.xlsx",sheet=2)

p1<-ggplot(data%>%filter(deadyear<2021),aes(x=deadyear,y=n))+
  geom_point()+
  geom_vline(aes(xintercept=2004), colour="red", linetype="dashed")+
  xlab("Year")+ylab("Number of deaths")+
  annotate("text", x = 2013, y = 48, label = "Year = 2004")+theme_bw()

p2<-ggplot(data%>%filter(deadyear<2021),aes(x=deadyear,y=alive))+
  geom_point()+
  geom_vline(aes(xintercept=2004), colour="red", linetype="dashed")+
  xlab("Year")+ylab("Number of living generals")+
  annotate("text", x = 2013, y = 1100, label = "Year = 2004")+theme_bw()

pdf("figure3.pdf", width=10, height=4)
ggpubr::ggarrange(p1,p2, ncol=2, nrow=1)
dev.off() 


####Table 2#######################################################

rm(list=ls(all=TRUE))

data<-read.dta13("1data_replicate.dta")
data$general_living_dummy<-0
data$general_living_dummy[data$general_living>0]<-1

model1 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow")
model2 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living+
                  lcp+loggdp+logpop+logarea+logdis_province+train2004+logdis_airport+ruggedness+
                  leader_incumbent+connect_prefecture_county+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow")
model3 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living_dummy+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow")
model4 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living_dummy+
                  lcp+loggdp+logpop+logarea+logdis_province+train2004+logdis_airport+ruggedness+
                  leader_incumbent+connect_prefecture_county+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow")
model5 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow",weights = inp)
model6 <- coxph(formula=Surv(`_t0`, `_t`, `_d`)~general_living+
                  lcp+loggdp+logpop+logarea+logdis_province+train2004+logdis_airport+ruggedness+
                  leader_incumbent+connect_prefecture_county+
                  cluster(countyid)+strata(prefectureid),data=data %>% filter(general_all>0),
                robust=T,method="breslow",weights = inp)

stargazer::stargazer(model1, model2, model3, model4, model5, model6, type = "text",
          keep = 'general_living',star.cutoffs = c(0.1, 0.05, 0.01),
          dep.var.caption  = "Time to 1st HSR Station",dep.var.labels.include = FALSE,
          keep.stat=c('n'),no.space=T,
          covariate.labels = c("No. of living generals", "No. of living generals (dummy)"))

#correct standerd errors
se = list(sqrt(model1[["var"]]), sqrt(model2[["var"]][1]),
          sqrt(model3[["var"]]), sqrt(model4[["var"]][1]),
          sqrt(model5[["var"]]), sqrt(model6[["var"]][1]))

####Figure A.3#######################################################

rm(list=ls(all=TRUE))
data<-read.dta13("1data_replicate.dta")

data<-data%>%filter(station==1)
graph<-as.data.frame(table(data$year),stringsAsFactors=F)
graph$Var1<-as.numeric(graph$Var1)

ggplot(graph,aes(x=Var1,y=Freq))+
  geom_bar(stat="identity",color="black", fill="white",alpha=0.1)+
  geom_text(aes(label=Freq), vjust=-0.3, size=3.5)+
  scale_x_continuous(breaks=seq(2005,2015,1))+
  xlab("Year")+ylab("No. of counties approved for construction")+theme_bw()
ggsave("figure.A3.pdf",width=8,height = 6)


####Figure A.4##############################################
#Figure A.4 is produced in QGIS


####Figure A.5##############################################

rm(list=ls(all=TRUE))

data<-read.dta13("1data_replicate.dta")

data<-data %>% filter(year==2004)

data$general_living_dummy<-0
data$general_living_dummy[data$general_living>0]<-1
table(data$general_living_dummy)

data$general_living_dummy <- factor(data$general_living_dummy, levels = c("1", "0"),
                                    labels = c("Living generals counties", "Non living generals counties"))

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(lcp))
result<-t.test(lcp ~ general_living_dummy, data = data)

p1<-ggplot(data, aes(x=lcp, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=0.7, y=18, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_x_continuous(breaks=seq(0,1))+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Least cost network counties",x="", y = "Density")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")


mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(loggdp))
result<-t.test(loggdp ~ general_living_dummy, data = data)

p2<-ggplot(data, aes(x=loggdp, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=4.5, y=1.2, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="County GDP",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(logpop))
result<-t.test(logpop ~ general_living_dummy, data = data)

p3<-ggplot(data, aes(x=logpop, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=1, y=1.5, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Population",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(logarea))
result<-t.test(logarea ~ general_living_dummy, data = data)

p4<-ggplot(data, aes(x=logarea, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=4.3, y=1.7, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Administrative area",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(logdis_province))
result<-t.test(logdis_province ~ general_living_dummy, data = data)

p5<-ggplot(data, aes(x=logdis_province, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=1.5, y=1.5, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3)))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Distance to province capital",x="", y = "Density")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")


mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(train2004))
result<-t.test(train2004 ~ general_living_dummy, data = data)

p6<-ggplot(data, aes(x=train2004, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=0.7, y=15, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"*"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_x_continuous(breaks=seq(0,1))+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Non-HSR station before 2004",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(logdis_airport))
result<-t.test(logdis_airport ~ general_living_dummy, data = data)

p7<-ggplot(data, aes(x=logdis_airport, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=4.5, y=1.8, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"*"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Distance to nearest airport",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(ruggedness))
result<-t.test(ruggedness ~ general_living_dummy, data = data)

p8<-ggplot(data, aes(x=ruggedness, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=3, y=1.4, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Geographical ruggedness",x="", y = "")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(leader_incumbent))
result<-t.test(leader_incumbent ~ general_living_dummy, data = data)

p9<-ggplot(data, aes(x=leader_incumbent, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=0.35, y=23, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3)))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_x_continuous(breaks=seq(0,1))+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Hometowns of incumbent leaders",x="", y = "Density")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(connect_pro_prefecture))
result<-t.test(connect_pro_prefecture ~ general_living_dummy, data = data)

p10<-ggplot(data, aes(x=connect_pro_prefecture, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=0.75, y=13, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3)))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_x_continuous(breaks=seq(0,1))+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="Provincial–prefecture party secretaries ties",x="", y = "Density")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(connect_prefecture_county))
result<-t.test(connect_prefecture_county ~ general_living_dummy, data = data)

p11<-ggplot(data, aes(x=connect_prefecture_county, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=0.7, y=14.5, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3)))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_x_continuous(breaks=seq(0,1))+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title=" Prefecture–county party secretaries ties",x="", y = "Density")+
  theme_bw()+theme(text=element_text(family="Times"),
                   legend.position="none")

get_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mu <- ddply(data, "general_living_dummy", summarise, grp.mean=mean(loggdp))
result<-t.test(loggdp ~ general_living_dummy, data = data)

p<-ggplot(data, aes(x=loggdp, color=general_living_dummy, fill=general_living_dummy)) +
  geom_histogram(aes(y=..density..), position="identity", alpha=0.2)+
  annotate(geom="text", x=4.5, y=1.2, 
           label=paste0("Mean Diff = ",round(result[["estimate"]][2]-result[["estimate"]][1],3),"***"))+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=general_living_dummy),
             linetype="dashed")+
  scale_color_brewer(name = "",palette = "Set1")+
  scale_fill_brewer(name = "",palette = "Set1")+
  labs(title="County GDP",x="", y = "")+
  theme_bw()+theme(legend.text=element_text(size=14),text=element_text(family="Times"),
                   legend.position="right")

p_legend <- get_legend(p)

pdf("figure.A5.pdf", width=16, height=10)
ggpubr::ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11, ncol=4, nrow=3, p_legend)
dev.off() 


####Figure A.6#####################################

rm(list=ls(all=TRUE))

county<-read_sf("./map/county.shp")
province<-read_sf("./map/province.shp")
line<-read_sf("./map/ldp.shp")
China <- st_read("./map/China.json")

data<-read.xlsx("2data_figure.xlsx",sheet=3)

data<-left_join(county,data)

data$leastdis[is.na(data$leastdis)]<-0
table(data$leastdis)

colourPalette <- c("#ffffff","#74c476")

ChinaMap<-ggplot() +
  geom_sf(data = China, colour = "black",fill=NA,lwd =0.5)+
  geom_sf(data = data, aes(fill=factor(leastdis)),colour = "#525252",alpha=0.8)+
  geom_sf(data = province,fill=NA,colour = "black",lwd =0.5)+
  geom_sf(data = line,colour="black",linewidth=1)+
  geom_sf(data = line,colour="white",linetype = "dashed",linewidth=0.5)+
  coord_sf(ylim = c(-2387082,1654989),xlim = c(-2387082,2554989),crs="+proj=laea +lat_0=40 +lon_0=104")+
  scale_fill_manual(name="",values=colourPalette,na.value = "white")+
  theme_minimal()+
  theme(legend.position="none",,
        panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

SouthChinaSea <- ggplot() +
  geom_sf(data = China, colour = "black",fill=NA,lwd =0.5)+
  coord_sf(xlim = c(117131.4,2115095),ylim = c(-4028017,-1877844),
           crs="+proj=laea +lat_0=40 +lon_0=104")+
  theme(aspect.ratio = 1.25,axis.text = element_blank(),axis.ticks = element_blank(),
        axis.title = element_blank(),panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour="black",fill = NA,linewidth=1.5),
        plot.margin=unit(c(0,0,0,0),"mm"))

ppi <- 600
png("figure.A6.png", width=10*ppi, height=6*ppi, res=ppi)
ggdraw() +
  draw_plot(ChinaMap) +
  draw_plot(SouthChinaSea, x = 0.8, y = 0.00, width = 0.19, height = 0.35)
dev.off()


####Figure A.7########################
rm(list=ls(all=TRUE))

coef1<-read.dta13("Coefficients1.dta")
coef1<-cbind(coef1,c(1:10))
names(coef1)[3]<-c("general_living")

pd <- position_dodge(0.8)

pdf("figure.A7.pdf", width=10, height=5)
ggplot(coef1 %>% filter(general_living<=7),aes(x=general_living,y=b))+
  geom_point(position=pd)+
  geom_errorbar(aes(min = b-1.96*se, max = b+1.96*se),position=pd,width=0)+
  geom_hline(yintercept =0,linetype=2)+
  xlab("No. of living generals")+ylab("Coefficient")+
  scale_x_continuous(breaks=seq(1,10,1))+
  scale_y_continuous(limits=c(-0.27,3.8),breaks=seq(0,3.5,0.5),labels=c("0","0.5","1","1.5","2","2.5","3","3.5"))+
  theme_bw() + theme(legend.position=c(0.85,0.7),
                     legend.background=element_blank(),
                     legend.key=element_blank(),
                     panel.grid.minor = element_blank(),
                     legend.text = element_text(size=14),
                     axis.text.x = element_text(size=14),
                     axis.text.y = element_text(size=14),
                     axis.title.x=element_text(face="bold",vjust=-.5,size=14),
                     axis.title.y=element_text(face="bold",vjust=1,size=14))
dev.off() 
